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Different approaches are presented to investigate diffusion from a point source in a slab delimited 
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consider the effect of absorption at the boundaries as well as the possibility that the particles that 
diffuse react with the diffusive medium. 
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I. INTRODUCTION 

The problem of steady-state diffusion from a source q{r) is described by a quite simple equation 

KAN{r) = -q{r) (1) 

where K is the diffusion coefficient (homogeneous to a surface per unit time). For a point source q{r) = S^(r) 
embedded in infinite space, the solution is also simple, given by 

For more complicated geometries of the diffusive volume, the solutions lose this simplicity. In some cases, it is helpful 
to use approaches different from brute-force resolution. We want to illustrate this point for the problem of diffusion in 
a slab delimited in the z direction by two absorbing boundaries, located at z = ±L and imposing that N(z = ±L) = 0. 
We present four approaches, giving identical results when the conditions of validity overlap, taking also into account 
the possibility that the particles that diffuse can also be destroyed by reacting with the diffusive medium. The equation 
then reads 

KAN - Tir}N = -q(r) (3) 

where the destruction rate r{r) is related to the density n{r) of the reacting medium, the reaction cross-section cr and 
the velocity v of the diffusing particle by T{r) = n{r)av. The authors first encountered this situation when studying 
the diffusion of cosmic rays emitted from sources in the galactic planoi. The galactic magnetic field has a stochastic 
component which is responsible for their diffusion, but also for their confinement in a volume which corresponds ap- 
proximately to the geometry described above. The destruction term corresponds to the nuclear reactions (spallations) 
that may occur when these cosmic Rays cross the regions of the galactic disk where the nuclei of interstellar matter 
are present. This is why we pay a particular attention to the situation where the sources, the measurement and the 
destruction process are localized in the plane z = 0. 

II. STEADY-STATE SOLUTION USING FOURIER-BESSEL TRANSFORMS 

A. Diffusion equation 

Given the geometry of the diffusion volume and the source, it is easier to use cylindrical coordinates. In all the 
following, the destruction term will not depend on the radial coordinate r, and we have 

(fiN IdN (fN T(z) q(r,z) 

! ! —N = - I (4) 

dr^ r dr dz"^ K K 
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FIG. 1: Geometry of the diffusive slab. The central plane contains the matter on which the diffusing particles can react. The 
upper and lower planes are absorbing boundaries, imposing a null density. 



We develop N{r, z) et q{r, z) over the Bessel functions Jolkr), i.e. 



N{r,z)= dkkJo{kr)N{k,z) 
'o 



N{k,z)^ / drr Joikr) N{r,z) 
Jo 



(5) 



(6) 



The radial dependence is now encoded in the relative weight of the Bessel functions. High values of k correspond to 
finer details in the radial distribution, much like for usual Fourier transforms. Inserting Q into l@J, and using the 
fundamental property of Bessel functions 



(7) 



we obtain 



dkk^Jo(kr)N{k,z)+ [ dkkJo(kr) \^^-^^^^^^-^^N(k,z)] =-— f dk k Mkr) q(k, z) (8) 
Jo \ dz^ K K Jo 



Using the property of orthonormalization 



c?rr Jo(fcir)Jo(fc2r) — 6(ki — fcj) 



(9) 



we select the equation for each mode k 

d^N{k,z) fT{z) 



dz^ 



K 



+ k^ N{k,z) 



qjk, z) 
K 



(10) 



B. Solution for a destructive plane 



If destruction is confined to the plane z = 0, we have T{z) = V5{z). We also consider a point-like source located in 
the disk, i.e. q{r, z) = S{z)6{Trr'^). The equation to be solved is, for each mode k, 



K 



K 



(11) 



In this expression, the combination T /K = l/?'d has the dimension of an inverse length, and the Bessel transform of 
a point-like source reads 



q[k) = I dr r Jo{kr) S{TTr ) = — 
Jo 277 



(12) 
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Outside of the disk {z ^ 0), the equation simplifies into 

The solution has to be even in z and satisfy the boundary conditions N{k, z = ±L) = 0. Therefore, it is given by 

The integration constant N(){k) is fixed by examining the solution in the plane z = 0. The derivatives of the solution 
(I14|l for z = are defined only in terms of distributions. The computation may be made easier using the identity 

exp |z| = exp(— z) + 20(z) sinh z (15) 

in the hyperbolic functions, as the derivative of the Heaviside distribution 0(z) is the Dirac distribution 5{z). This 
yields 

N{k, z ^ 0) = iVo(fc) _ 2No(k)0{~z) coth(fcL) sinh(yfcz) (16) 

smh {kL} 

The second derivative reads 

d^N{k,z) 
d? 

Inserting this last expression into Hll|) yields 



fc2 N{k, z) - 2Na{k)k5{z) coth(A;L) (17) 



and the final solution is 



kMkr)dk sinh{fc(L-|z|)} 
' 7o 2TT{T + 2Kkcoth{kL)) sinh{fcL} ^ ' 

The diffusion process acts as a filtering in the Bessel space, the diffused density being related to the source by the 
transfer function 

^^^^ " T + 2KkcotHkL) 

It is a low-pass filter. Diffusion tends to erase the small scale features. The solution in the disk (z — 0) can also be 
written under a form making apparent the correction to the free diffusion case ^ 

^' ' 4TrKrJ„ r/2rd + fcr coth(fcL) 

This is illustrated in Fig. [3 



C. Alternative formulations for a better convergence 

The integral involved in Eq. H21() is of the form 

/>OC 

/[/] = / Mx)f{x)dx (22) 
Jo 

When /(oo) = 1, as is the case up to a normalization in Eq. H21() when z = 0, the slow decrease of the oscillations in 
the integrand makes the numerical computation of /[/] quite tricky. A few examples of manipulations which help to 
improve the numerical implementation of Eq. (|21fl are given in the appendix. 
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FIG. 2; Deviation from the free diffusion case l/47ri('r, as a function of r/L for the case without reaction (thick curve: finite 
L, r = 0) and of r/r^ for the case without boundary (thin curve: L — > oo, finite F). 



D. Solution for a short-lived species 

The same equation is also relevant to study the diffusion of unstable particles, with a lifetime r. The destruction 
rate is then homogeneous inside the diffusive volume, and given by = 1/r. It is straightforward to show, using the 
same procedure as above, that in that case 



M( 1 r krJoikr)d{kr) /.2 ^ 



III. TIME-DEPENDENT SOLUTION 



A more precise description of the diffusion process may be obtained from the time-dependent diffusion equation. 
This gives another formulation of the steady-state solution, as a series having better convergence properties. It also 
enables to take differently into account the case of decaying particles. 

A. The time-dependent diffusion equation 

The time-dependent diffusion equation reads 

ON {I d f dN\ d^N] 

Wc seek the solution at t > such that N{r,z,t = 0) = S^{r). It is convenient to use the typical length = K/T 
introduced above, or alternately its reciprocal fed, so that Eg. 1241 mav be written as 

dN {I d f dN\ a^iVl , 

+^ -^d^W^ (25) 



d{Kt) \r dr \ dr J dz^ 
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The diffusion process in the z direction and in the radial direction arc independent. As only pure diffusion occurs in 
the radial direction, the density can be written as 

AA(r,^,t) = ^e-^'/^^*iV(z,i) (26) 

where the function N{z,t) satisfies the time dependent diffusion equation along z 

dN d^N 



d{Kt) (9z2 



- k^5{z)N (27) 



B. Derivation 

First, we seek solutions of the form N{z, t) = f{z)g{t), which separates the diffusion equation into 

g' = ^ag and - af = Kf" -TSiz)f (28) 

where a must be positive in order to eliminate runaway solutions for g. The equation on / can be solved for z > 
and z < with the condition f{±L) = as 

f{z) = Asm{k{L-\z\)} (29) 



where k = y^a/K. Derivation in the sense of the distribution, as in the previous section, yields 

Kf"{z) = -af{z) - 2KkAS{z) cos kL (30) 

Inserting into H28|l gives the condition 

2A: cotan(A;L) = — fcj (31) 

There is a infinite discrete set of A:„ satisfying the above condition, which gives the allowed values a„ = Kk^ of a. 
For example, whit no reaction (kd = 0), kn = {2n + 1)tt/2L. The general solution reads 

oo 

N{z,t) = ^ A„e-""*sin{fc„(L- |z|)} (32) 

71=1 

The functions sin{fc„(L — \z\)} form an orthogonal set, and it is found that 

nL 

sin{/c„(L - \z\)}sin{kn'{L - |z|)} dz = (5„„'C„ (33) 



with 



sin2fc„L 

Cn = L- —— (34) 



The An are found by imposing that for t = 0, the distribution is a dirac function, 

oo 

5(z) = ^A„sin{fc„(i- |z|)} (35) 

n=l 

Multiplying by sin{fc„i(L — |z|)} and integrating over z yields 

Am = sin kjnL (36) 

so that finally 

oo 

N{z, t)^Y^ c-^e"""* sin(fc„L) sin {fc„(L - |z|)} (37) 

and 

= i^^^'P (-^^ j E^"'^"""*^i^(^»^) sin{/c„(L- |z|)} (38) 



n=l 



The radial distribution in the disk is given by 



^ = 0' = rio {^iw; ) E ^.T'e-^-^* ^in^k^L) (39) 



2 \ °° 



Tl=l 



FIG. 3: Particle distribution as a function of zjL and r/L, for Kt = 0.05 L and Kt = 0.2 L. At early times, the distribution is 
close to the free case, as very few particles had time to reach the boundary. At later times, the effect of absorption are more 
pronounced. 



C. Reformulation of the steady-state model 

The stationary regime results from the contimious superposition of solutions for instantaneous sources, so that the 
corresponding solution is given by 



A4tat(r,z)= / J\f{r,z,t)dt (40) 
Jo 

Using the identity (Gradshteyn & Ryzhik 1980) 

^ e-«t-/3/t = 2Ko (2^^) (41) 
involving the Bessel function of the third kind Kq, the integration of H38|) yields 

^ oc 

Mtat(r,z) = 7r-ry^c~^Ko (fc„r) sin(fc„L) sin{fc„(L- z)} (42) 

n=l 

where the Bessel function of the third kind Kq has been introduced. The density in the disk is thus given by 

Mtat(r, z = 0) = -—Y, c-'Ko (fc„r) sin^ik^L) (43) 



2'kK 

n=l 

This expression provides an alternative (but is exactly equivalent) to the usual Fourier Bessel expansion, using the Jq 
functions. It is particularly well suited for sources well localized in space, like point-like sources, because the functions 
over which the development is performed (the Kq) do not oscillate. As a consequence, convergence of the series above 
is fast and the expression above provides a powerful alternative to compute the density for not too small values of 
/L, as illustrated in Fig. Q for the case F = 0, for which sin^ fc„L = 1 and c„ = L, so that 



r 



^■„,,..„,._L_fK„f?2±ii;i) (44) 



2TrKL ^ " V 2 L 
This case is illustrated in Fig. Q). 

D. The case of unstable particles 

When the particles can decay with a rate F^ = l/r, the previous time-dependent expressions must simply be 
corrected by a multiplicative factor exp(— F„t), which amounts to make the substitution — > Q!„ -t-F^ so that finally, 
the expression H43|) is still valid provided the substitution fc„ — > \/kf-, +Tu/ K is performed. 
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FIG. 4: Values of the correction to the free diffusion case as a function of r/L, for F = 0, computed with the series in Kq, 
truncated at different numbers of terms. The first term alone reproduces quite well the profile for r > 2. 



IV. METHOD OF IMAGES 



A. Definition and result 



In the absence of destruction in the disk, a completely different approach is provided by the elegant method of the 
images inspired from electrostatics. In infinite space, the solution to our diffusion problem is straightforward and is 
given by the relation Q in the case of an initial Dirac distribution at the source S. We furthermore would like to 
impose the boundary condition according to which the cosmic-ray density N vanishes on the plane z = +L at any 
time. To do so, we can introduce the virtual source S" that is the image of the real source S with respect to the 
boundary z = +L acting like a anti-mirror. The cosmic-ray densities which S and its image S' generate are equal up 
to a relative minus sign that allows both contributions to cancel out exactly on the boundary. In order to impose that 
the density also vanishes at z — — L, we can consider the anti-image S" of S with respect to that lower boundary. 
Because two anti-mirrors are now present at z = +L and z ~ — L, an infinite series of multiple images {Sn} arises. 
They are aligned with the real source Sq = S in the vertical direction and the position of Sn is given by 

s„ = 2L?i+ (-l)"s . (45) 

The virtual source S'„ results from \n\ reflections throughout the mirrors and its production is affected by a sign (—1)" 
with respect to the real source S. The distribution of sources within the Galaxy is not perturbed by the presence of 
their virtual images that are located outside the domain of interest. 

We readily infer that in the presence of boundaries, the time-dependent solution 10) is modified into 



n— — oo ^ ^ 

which can be rewritten, decoupling the diffusion processes along the vertical axis and in the radial direction, as 

„-r^/4:Kt ( +°° (T\n , 1 



71— — OO 
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FIG. 5: Position of the virtual sources, images of the real sources through the anti-mirrors nt z — ±L. The black dots represent 
positive sources (they contribute like the real source) and white dots represent negative sources. 



The steady-state solutions given by 



N,{r,z,t)= J2 , /, \ (48) 



11— — oo 




FIG. 6: Density in the plane z = 0, as a function of r/L, and normalized to the free diffusion case. The dots represent the 
exact solution, and the thin lines show the result given by the series of images, for several truncatures of the sum. At very 
large distances (r ^ L) and in the case of a small number of sources, these are felt as a single source (positive or negative, 
depending on the parity of TV), so that the density tends to plus or minus the free diffusion case. 



B. Equivalence with the first approach 



The expression (|21ll can be easily transformed into (|48|1 . for z = in the case T = 0. Indeed, in the expression 

1 

N(r,z)^ / Jo(fcr) tanh(fcL)(i(fcr) (49) 

AnKr Jq 

using the development 

oc 

tanh(fcL) = (1 - (l + e'^*^'^) = (l - g-^'^i) ^(_1)" e'^^'^^ (50) 
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a few steps lead to 



-i- 1 + 2^:/ Jo( 



N(.r, z) ^ ^1 + 21^/ Mkr) (-1)" e^^^^ d(fcr) ^ (51) 



Which finally gives the expression H48() . using Jo{ax)e exp{—bx)dx = (a^ + 5^) 

C. Equivalence with the former approach 

The diffusion from the single source S within a slab - on the boundaries of which the cosmic-ray density vanishes - 
amounts simply to the diffusion in infinite space from the series of sources S'„. Along the vertical direction, the initial 
distribution is 

7V,(z,0= (-ir^(^-s„) . (52) 

n— — oo 

Because that distribution is periodic, its Fourier transform which we define by 

/ + 00 
n{s,0->z,0} e-'^^'dz , (53) 
-oo 

is composed of a discrete series of modes k. It may actually be expressed as 

+ 00 ^ +00 ~j 



n— — OO I n— — oo 



After some straightforward algebra, very similar to what is done in optics to compute the diffraction from an infinite 
diff"racting grid, the sum S is transformed into 

+ 00 

n— — oo 

We conclude that the initial vertical distribution n{s,0 z,0} contains a series of modes with discrete wavevectors 
k = n7r/2L. Inserting this expression into H54|l . we see that odd values of n ^ 2p — 1 are associated to modes for 
which kpL = {p — 1/2) tt and contribute a factor 

7r 

N {kp) = — cos kpS , (56) 

whereas even values of n = 2p lead to wavevectors k'p such that k'pL = pi: and to 

N {k'p) ^ sill k'ps . (57) 

The initial cosmic-ray distribution in the vertical direction n {s, — > z, 0} may be expressed as a Fourier series on the 
various odd and even modes p which we have just scrutinized 

r I 1 1 

Ns{z, t) = ^ < — cos {kps) cos {kpz) + — sin {k'ps) sin (fc^z) > . (58) 
p=i ^ ^ 

Because each Fourier mode k exponentially decays in time like exp [—Kk^t) as a result of diffusion, the initial 
distribution subsequently evolves into 

^(-Kk^t ^-Kk'lt ] 

Ns{z, t) = N < cos {kps) cos {kpz) H sin {k' s) sin (fc' z) > . (59) 

P^il ^ ^ J 

When the source is iat the origin (s = 0), the density at the origin (z = 0) is given by 

+00 ^-Kkpt 

n{s,0^z,t} = ^ — cos (fcpz) . (60) 

p=i 

which, when radial diffusion is taken into account, is equivalent to the expression H38|l in the case of no reaction (then 
sinfc„i = 1 and cosfc„L = 0). 



10 




z(t)>L 







FIG. 7: For each path such that max > L, there are two paths, one of which satisfies z{t) > L but not the other. 

V. RANDOM WALK APPROACH 

It is well known that diffusion is closely related to random walks. The density at each point is related to the number 
of stochastic path which reach this point, from the source. When a boundary is present at z = ±L, the paths that 
would go beyond this boundary do not contribute anymore to the density and must be discarded. When destruction 
can occur in the disk, the paths that cross the disk should be attributed a lower weight. In this section, we investigate 
separately these two effects. 



A. Probability of not escaping 

We first compute the probability that a particle emitted at time t = in the disk has not reached the boundary at 
a further time t. At this time, position z is given by a random walk of duration t. The probability we seek is given by 

{(max < L) n (min > -L)} = 1 - "P {(max > i) U (min < -i)} (61) 

where max and min are the maximum and minimum z reached by the random walk. From elementary statistics, 

V {(max > L)U (min < -L)} = V {max > Lj+V {min < -L} - T' {(max > i) n (min < -L)} (62) 

By symmetry, V {max > L} = V {min < —L}. Furthermore, using the principle of reflection, it is straightforward to 
show that (see figure [71) 



V {max > L} = 21^ {z{t) > L} 

It is also apparent that (see figure (SJ 

P {(max > L) n (min < -L)} = V {(max > 3L) U (min < -3i)} 

so that finally 

P {(max < L) n (min > -L)} = I - AV {z(t) >L}+'P {(max > SL) U (min < -3L)} 
Now, the same reasoning with Eq. H62() applied to the last term yields 

V {(max <L)r\ (min > -L)} = I - AV {z{t) > L] + AV {z{t) > 3L} - AV {z{t) > 5L} + 



(63) 
(64) 
(65) 
(66) 



B. Probability of not escaping for particles reaching the disk 

The problem we addressed in the previous sections was a bit different, though, as we were interested in particles 
reaching a given point, for example in the disk. The same reasoning as before can be used, and the probability that 
a particle reaching the disk at time t has not wandered farther than the boundaries is given by 



Vd {(max < L) n (min > -L)} = 1 - Vd{{ma.x > i) U (min < -L)} 



(67) 



3L - 







max >L & min <-L 



FIG. 8: For each path such that max > L and min < — L, there is a symmetric path satisfying z(i) > 3L. 
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FIG. 9; If a path such that z{t) — wanders beyond the upper boundary, then there is a symmetric path that satisfies z{t) 
The probabihty that the path goes beyond this boundary is thus given by the ratio rapport V {z{t) — 2L} jV {zii) — 0} 

where we have used the compact notation for conditional probabihty 

Vd {event} = V {event |z(t) = 0} 
Exactly as before, the use of symmetry yields 

Vd {(max > L)U (min < -L)} = 2Vd {max > L} - "Pd {(max > L) n (min < -L)} 
Now, the principle of reflection yields a different result for the conditional probability (see figure 

r{z{t) = 2L} 



As before. 



so that finally 



7'{max>L|.(0 = 0}= 7'{.(t).0} 
Vd {(max > L)n (min < -L)} = Vd {(max > 3L) U (min < -3L)} 



Vd {(max < L) n (min > -L)} = 1 - + T'd {(max > 3L) U (min < -3L)} 

7'{z(i) = 0} 
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The same reasoning, applying equation H69|l to the last term, gives 

X . M .X 1 V{z(t)=2L} V{z(t) = 4L} r{z{t) = 6L} , , 

{(max < L) n (mm > -L)W) ^ 0} = 1 - 2^LU^ + ^^Who/ " ^ V {zit) = 0} + ' ' ' ^^^^ 

or 

P{(max < L)n (min > -L)\z{t) = 0} = 1 + 2^(-l)" exp ( ) (74) 



It has the same form as given by the method of images. 



C. The effect of destruction in the disk 



We now investigate the effect of destruction occurring in the disk, in the simple case where no boundary is present. 
As before, the steady-state 3-D solution will be obtained from the one-dimensional time-dependent diffusion. We 
denote p the probability that a particle crossing the disk is destroyed. In the case of a random walk z = X^t^i '^i 
consisting of N = t/r elementary steps of length A and duration t, the probability distribution of disk-crossing 
numbers n is given by (Papoulis 2002) 

2.TIT ( TI^t\ 

drdin\t) = dV{n\t, z{t) = 0) = exp dn . (75) 

K2t \ K2t J 

In this expression, K2 ^ 1 depends on its statistical properties (for instance, K2 = 2 for elementary steps Zi = ±A and 
K2 ~ 1.43 for Zi uniformly distributed in the interval [—A, A]). The diffusion coefficient is defined as 

^-^(f!)^iV((^A^_^ (76) 
2t 2N T ^ 2 ' ^ ' 

where K3 = {{zij X)^) is the variance of the elementary random step (in units of A) and v = A/r. We thus finally have, 

,^ / , ^ '^Kn ( 2Kr? \ , 

dVd{n\^) = 5- exp dn (77) 



K2K3VH \ K2K^V^t 

We are now able to compute the probability distribution of disk crossings for particles emitted from a distance r in 
the disk as 

r'^Vdmdt. (78) 
dn Jq dn 

where the probability that a CR reaching distance r in the disk was emitted at time t is 

^^WO-(]^exp(-^) . (79) 



The above integral (|78() can be performed, yielding the final result 

'1 + # ) , (80) 



dVd{n\r) rln ^ .2.2 x -3/2 



dn 

with 7'q = 8/C2/k2K3i;2 ~ 2\^k,^/k2. We can also compute the integrated probability, that more that riQ crossings 
have occurred, as 

Va{n>n^\r)=(^+'^^ ' . (81) 

A particle having crossed n times the disk has the probability p„ = (1 — p)" ^ exp(— np) of surviving, so that the 
survival probability at distance r is given by 

Jo dn 

dx ^-xrp/ro 



{l + x-^f^ 
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The density of Cosmic Rays in the disk is then given by 



/{r) ^1 I" X dx ^-xrp/ro 



N(r) = = / g-xrp/ro (82) 

AnKr AnKr Jq [I + x^) ' 
This form is actually equivalent to H19(l , with L ^ oo, 

Indeed, rewriting l/27r(r + 2Kk) as dy cxp(— 27r(r + 2Kk)y) and reversing the integrations, (|21fl can be written, 

poo pOO 

N{r,z = 0)= dye-^''^y k Jo{kr) dke-'^'^^'^y (84) 



Using the identity (Gradshteyn & Ryzhik 1980) 

(q,2 + ^2)3/2 



e-^yjoiPy) y dy = , ^ , (85) 



we finally have, performing the change of variables AirKy / 



r = X, 



OO 



^ ' ^ io ^ (r2 + (47rX2/)2)3/2 ^ AirKr {l^x^f/^ ^ ' 

This equation has the same form as H82|l . and relates the microscopic and macroscopic properties of diffusion, as p/rQ 
must be equal to T/ttK. It should also be remarked that this integral is easier to compute than those involving Jq 
functions, having a faster convergence as the integrand does not oscillate. 



VI. DISCUSSION 



The four methods presented in this paper do not all apply in every situation. The first two are valid for arbitrary L 
and r as well as with spontaneous decay of unstable species. The second one contains a richer physical information, as 
it provides the density as a function of time. The third one is only valid when F = 0, taking only into account the effect 
of absorption by the boundaries. The fourth one is more general as we have also presented the case F ^ 0. However, 
it needs to be worked further out in order to consider simultaneously absorption at the boundaries and destruction 
in the central planed. Comparison of the method of images and the use of random walk give some complementary 
insights to the consequences of absorption at a boundary. The paths that wander beyond the boundaries, and that 
we had to suppress by hand in Sec. are actually those connecting the point at which the density is sought, to the 
negative images introduced in Sec. IIVI The effect of these negative images is to destroy the paths that would wander 
out of the diffusive volume. It is important to note that this is different from the path integral interpretation of the 
Schrodinger equation, in which case the paths are weighted by a complex phase term, whereas in the classical diffusion 
case we have discussed, the paths have only ±1 factors. 
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APPENDIX A: TRANSFORMATIONS OF EQ. (21) FOR AN EASIER NUMERICAL COMPUTATION 

a. Substraction Using 

/ Jo{x)dx = 1 (Al) 
"'0 

/[/] may be rewritten 

/•OO 

I[f] = l- Mx)il- f{x))dx . (A2) 
Jo 

The convergence is faster, as 1 — / ^ when a; ^ 00. 
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b. Integration by part Using the identity (xJi)' = xJq and integrating by parts, one has 

/[/] = - hm [Ji(x)/(x)] + r Ji(x) f ffl - fix)) dx . (A3) 

x^0+ Jo \ X J 

Using the identity Jq = — Ji and integrating by parts again, 
/[/] = - hm [Ji(x)/(x)] + 

These expressions provide several efficient alternatives to evaluate /[/], provided the integrated terms are well defined, 
i.e. if 

f{x) = A + Bx\nx + 0{x) (A5) 

c. Comparison to a known function Part of the difficulty to evaluate numerically the Bessel expansions comes 
from the fact that the original functions are singular at the source position. As a result, the large k modes continue 
to be important to reconstruct the solution. We can take advantage of the fact that the singularity is known, as the 
density is quite close to the free diffusion case fici{r) = l/AirKr for small r. The corresponding Bessel coefficients are 
given by 

^^^W=fr^^Mkr)dr^^^ (A6) 

It is then judicious to write the density as 

N{r, 0) - /,ef(r) + Mkr) {/(fc) - /ref(fc)} dk = + Jo(fcr) jfc /(fc) - ^| rffc (A7) 

where the singularity is entirely contained in the first term: the Bessel expansion has been regularized. In the 
remaining integral to be computed, the convergence is much faster, as the large k modes contribute very little. Other 
choices of /rof('') may be preferred for particular values of L, T and r. This method yields a very good and rapid 
convergence for sources located in the thin disk z = 0. 

d. Softening of the source term Finally, the source term may be spread out on a radius a, by replacing the point 
source S{r) by a disk source q{r) = 0{a — r)/{'Ka^) for which an extra 2Ji{ka)/ka term appears in the Bessel transform. 
With a judicious choice of the parameter a, the solution is very close to the original for r 3> a, but convergence is 
much faster due to the extra 1/k factor. 



Mx)[^~nx) 



+ Mx) f"{x) 



dx 



(A4) 
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^ For example, the probability distribution of disk crossing can be found in Taillet et al., |astro-ph/0308141| for a general case 
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FIG. 10: Values of the correction to the free diffusion case, for r = 0.01 L and F = 0, computed with the different methods 
described in the text, for different values of the upper boundary in the integral. 



